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Abstract. We have carried out a multifrequency analy- 
sis of the radio variability of blazars, exploiting the data 
obtained during the extensive monitoring programs car- 
ried out at the University of Michigan Radio Astronomy 
Observatory (UMRAO, at 4.8, 8, and 14.5 GHz) and at 
the Metsahovi Radio Observatory (22 and 37 GHz). Two 
different techniques detect, in the Metsahovi light curves, 
evidences of periodicity at both frequencies for 5 sources 
(0224 + 671, 0945 + 408, 1226 + 023, 2200 + 420, and 
2251 + 158). For the last three sources consistent peri- 
ods are found also at the three UMRAO frequencies and 
the Scargle (1982) method yields an extremely low false- 
alarm probability. On the other hand, the 22 and 37 GHz 
periodicities of 0224 + 671 and 0945 + 408 (which were 
less extensively monitored at Metsahovi and for which 
we get a significant false-alarm probability) are not con- 
firmed by the UMRAO database, where some indications 
of ill-defined periods about a factor of two longer are re- 
trieved. We have also investigated the variability index, 
the structure function, and the distribution of intensity 
variations of the most extensively monitored sources. We 
find a statistically significant difference in the distribu- 
tion of the variability index for BL Lac objects compared 
to flat-spectrum radio quasars (FSRQs) , in the sense that 
the former objects are more variable. For both popula- 
tions the variability index steadily increases with increas- 
ing frequency. The distribution of intensity variations also 
broadens with increasing frequency, and approaches a log- 
normal shape at the highest frequencies. We find that vari- 
ability enhances by 20-30% the high frequency counts of 
extragalactic radio-sources at bright flux densities, such as 



Send offprint requests to: A. Ciaramella: ciaram@unisa.it 



those of the WMAP and Planck surveys. In all objects 
with detected periodicity we find evidence for the exis- 
tence of impulsive signals superimposed on the periodic 
component. 

Key words: Methods: data analysis - BL Lacertae ob- 
jects: general - quasars: general - radio continuum: general 



1. Introduction 

The name Blazars identifies a family of radio-loud Ac- 
tive Galactic Nuclei (AGNs) showing a rather complex 
phenomenology: extreme variability at all wavelengths, 
polarization, strong 7-ray emission, brightness tempera- 
tures exceeding the Compton limit ( |Urry 1999| ) . The large 
amount of work done in the last decade has led to a rather 
general consensus on the global mechanism responsible 
for the emission: a rotating black hole surrounded by a 
massive accretion disk with an intense plasma jet closely 
aligned to the line of sight. Relativistic electrons produce 
the soft photons through synchrotron emission, while hard 
photons are produced by inverse Compton scattering. This 
overall scenario, however, still presents a large number 
of poorly understood details which, in turn, lead to a 
wide variety of models and call for long term and multi- 
wavelength campaigns capable to provide the necessary 
observational constraints. 

Variability measurements provide key information on 
the AGN structure, down to linear scales or flux density 
levels not accessible with interferometric imaging. The 
most extensive blazar monitoring campaigns have been 
carried out at radio frequencies. The University of Michi- 
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gan monitoring program at 4.8, 8.0 and 14.5 GHz has ob- 
tained data on over 200 sources for over three decades. At 
higher frequencies, the Metsahovi group (Terasranta et al. 
1998) have reported observations at 22, 37 and 87 GHz of 
157 extragalactic radio sources, many of which have been 
monitored for over 20 years. The Bologna group (Bondi et 
al. 1996) have observed at 408 MHz 125 radio sources for 
15 years. A multifrequency monitoring of large sample of 
compact radio sources has been carried out by Kovalev et 
al. (2002) with the RATAN-600 telescope. 

The mechanisms for variability are still not well un- 
derstood. Possibilities discussed in the literature include 
shocks in jets (Marschcr & Gear 1985; Aller, Aller & 
Hughes 1985; Marscher 1996) and changes in the direc- 
tion of forward beaming [due, e.g. to helical trajectories 
of plasma elements (Camenzind & Krockcnberger 1992) 
or to a precessing binary black-hole system (Begelman et 
al. 1980; Sillanpaa et al. 1988)], introducing flares due to 
the lighthouse effect. Thus, variability furnishes important 
clues into size, structure, physics and dynamics of the ra- 
diating source region. 

The cm- wavelength variability of extragalactic sources 
has been investigated by Hughes et al. (1992) and Aller 
et al. (1996, 2003) exploiting the uniquely large, long 
term, University of Michigan Radio Astronomy Obser- 
vatory (UMRAO) data-base. Kelly et al. (2003) applied 
a cross-wavelet transform to analyze the UMRAO light 
curves of the Pearson- Readhead VLBI survey sources. The 
Metsahovi data-base has been analyzed by Lainela & Val- 
taoja (1993), Terasranta & Valtaoja (1994), Valtaoja et 
al. (1999), Lahteenmaki et al. (1999), Lahteenmaki & Val- 
taoja (1999). 

Still, the information content of the radio-monitoring 
data-base has not yet been fully exploited. We are inter- 
ested, in particular, in utilizing it to predict the effect of 
radio-source variability on the high-frequency all-sky sur- 
veys to be carried out by ESA's Planck mission. Such 
predictions are useful to plan the quick-look analysis of 
Planck data and to organize follow-up observations. To 
this end, the Metsahovi data, taken at frequencies close to 
those of Planck-LFI (Low Frequency Instrument) chan- 
nels, are particularly well fit . 

Another very interesting, highly debated issue, is the 
possible presence of periodicities in the blazar light curves. 
Claims for the existence of periodic behavior arc mostly 
based on optical data (Lainela et al. 1999; Fan 2000, Fan 
et al. 2002, and references therein). The most famous, but 
still somewhat controversial, case is the 11.65 year peri- 
odicity of OJ287 (Sillanpaa et al. 1988, 1996; Marchenko 
et al. 1996; Hagen-Thorn et al. 1997; Pietila, et al. 1999). 
Evidence of a persistent modulation of the radio total flux 
and polarization of this source, with period of ~ 1.66 yr 
was found by Aller et al. (1992) and Hughes et al. (1998). 

A possible periodicity of ~ 5.7 years in the radio light- 
curve of the BL Lac object AO 0235 + 16 has been reported 
by Roy et al. (2000) and Raiteri et al. (2001). Indications 



of a periodic or quasi-periodic component with a time- 
scale of 5.5-6 years (close to that of AO 0235 + 16) of the 
radio emission of the BL Lac object S5 0716 + 714 were 
found by Raiteri et al. (2003) 

In this paper we present a new investigation of ra- 
dio light curves of a sample of blazars, exploiting more 
effective techniques than the Pcriodogram method com- 
monly used (in various versions) for periodicity analysis. 
The maximum time interval covered by the data we used is 
1979-2001 in the case of Metsahovi, 1965-1999 in the case 
of UMRAO. The techniques used are described in Section 
2 and the main results are presented in Section 3. In Sec- 
tion 4 we discuss the multifrequency variability properties 
for a blazar sample monitored by both the UMRAO and 
the Metsahovi Radio Observatory, and estimate the effect 
of variability on high-frequency counts of such objects. In 
Section 5 we summarize our main conclusions. 



2. The time series analysis 

In what follows, we shall assume x to be a physical variable 
measured at discrete times ti\ x(U) can be written as the 
sum of the signal x s and random errors R: 

x l = x(U) = x s (U) + R(ti) . (1) 

The problem is how to estimate possible fundamental fre- 
quencies which may be present in the signal x s (ti). In the 
case of even sampling, many Fourier-like tools can be ef- 
fectively used (cf. Fernie 1979; Home & Baliunas 1986). 
These methods, however, encounter problems when deal- 
ing with unevenly sampled data and shortcuts like the 
use of interpolation to re-sample the data, introduce a 
noise amplification which usually undermines the subse- 
quent use of Fourier based techniques which are notori- 
ously very sensitive to the noise level of the input data 
(cf. Kay 1988; Marple 1987, but see also Horowitz 1974). 

2.1. The Periodogram (P) method 

The most commonly used tool for periodicity analysis of 
both evenly and unevenly sampled signals is the so called 
Periodogram (hereafter P), which is an estimator of the 
signal energy in the frequency domain (Deeming 1975), 
and has been extensively applied to the study of light 
curves of variable stars, both periodic and semi-periodic. 

In the case of even sampling, the Periodogram has a 
simple statistic distribution (exponentially distributed for 
Gaussian noise). This, however, is no longer true for un- 
evenly sampled data, thus making it difficult to control er- 
rors (Kay 1988; Marple 1987; Oppenheim & Shafler 1965). 
Scargle (1982) and Lomb (1976) introduced a modified 
form of Pcriodogram which takes explicitly into account 
the effects of uneven sampling. 

Let / be the frequency and r a shift variable, and let 
us also suppose that we are dealing with an observation 
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Table 1. Objects for which we have positive detection of periodicity on Metsahovi daily averaged data. Column 1: 
object identification. Columns 2-5: 22 GHz data; columns 6-9: 37 GHz data. Columns 2 & 6: number of data points; 
columns 3 & 7: maximum admissible period to avoid aliasing; columns 4 & 8: period (in units of 10 3 days) obtained 
by STIMA; column 5 & 9: period obtained from the Lomb's Periodogram. 



Name 


N 


Mx. P. 


STIMA 


Lomb 


N 


Mx. P. 


STIMA 


Lomb 


Notes 






(xlO 3 ) 


(xlO 3 ) 


(xlO 3 ) 




(xlO 3 ) 


(xlO 3 ) 


(xlO 3 ) 




(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


(9) 


(10) 


0224 + 671 


76 


3.062 


1.021 


1.021 


63 


3.668 


0.950 


0.970 




0945 + 408 


47 


3.742 


1.386 


1.336 


26 


4.467 


1.313 


1.240 


few points 


1226 + 023 


694 


6.665 


3.029 


2.777 


716 


7.822 


3.260 


1.261 




2200 + 420 


644 


6.676 


3.034 


3.034 


715 


7.873 


2.811 


3.028 




2251 + 158 


571 


6.676 


2.384 


2.384 


549 


7.538 


2.217 


2.512 





Table 2. Analysis of UMRAO data on Metsahovi sources with detected periodicity (Table 0). Columns 2-5, 6-9, and 
10-13 refer to 4.8, 8.0, and 14.5 GHz data, respectively. 



Name 


N 


Mx. P. 


STIMA 


Lomb 


N 


Mx. P. 


STIMA 


Lomb 


N 


Mx. P. 


STIMA 


Lomb 






(xlO 3 ) 


(xlO 3 ) 


(xlO 3 ) 




(xlO 3 ) 


(xlO 3 ) 


(xlO 3 ) 




(xlO 3 ) 


(xlO 3 ) 


(xlO 3 ) 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


(9) 


(10) 


(11) 


(12) 


(13) 


0224 + 671 


110 


5.998 


1.764 


1.764 


119 


6.859 


2.286 


2.286 


110 


7.592 


2.109 


2.109 


0945 + 408 


73 


6.048 






126 


7.993 






92 


5.920 


2.690 


2.466 


1226 + 023 


493 


6.424 


3.197 


3.197 


1294 


1.232 


3.082 


3.082 


760 


9.114 


3.038 


3.038 


2200 + 420 


745 


7.844 


2.801 


2.801 


1212 


1.131 


2.261 


1.413 


1056 


9.056 


2.830 


2.830 


2251 + 158 


554 


7.674 


2.240 


2.240 


1175 


1.195 


2.390 


2.490 


885 


9.190 


2.297 


2.297 



Table 3. Summary of detected periods (in units of 10 3 days). 



Name 


4.8 GHz 


8.0 GHz 


14.5 GHz 


22 GHz 


37 GHz 


0224+671 


1.764 


2.286 


2.109 


1.021 


0.950 


0945+408 






2.690 


1.386 


1.313 


1226+023 


3.197 


3.082 


3.038 


3.029 


3.260 


2200+420 


2.801 


2.261 


2.830 


3.034 


2.811 


2251+158 


2.240 


2.390 


2.297 


2.384 


2.217 



2.2. Autocorrelation matrix based analysis: the STIMA 
approach 

Recently, more effective techniques based on the analysis 
of the signal autocorrelation matrix were introduced [for 
a detailed exposition of the problem see Kay (1988) and 
Marple (1987)] . In a previous paper (Tagliaferri et al. 1999, 
hereafter T99) some of us introduced the so called STIMA 
algorithm, based on a particular type of the MUlty Signal 
Classifier (MUSIC) (Oppennheim & Schafer 1965) esti- 
mator specifically tailored to work with unevenly sampled 
data and on a robust nonlinear PCA Neural Network used 
to extract the principal components of the autocorrelation 
matrix of the input sources. Without entering into details 
(which may be found in T99) we now briefly summarize 
its main features. 

Let us assume to have a signal with p sinusoidal com- 
ponents (narrow band). The p sinusoids are modelled as 
a stationary ergodic signal, and this is possible only if the 
phases are assumed to be independent random variables 
uniformly distributed in the interval [0,27r]. To estimate 
the frequencies of these components we exploit the prop- 



series formed by N points x{n). Their mean and variance 
are given by: 

1 N 

* = N £ X{n) ^ = £»=i ~ ■ (2) 



The normalized Lomb's P L , i.e. the power spectrum as a 
function of the angular frequency lu = 2irf > is defined 
as: 



1 

2^2 

1 

2^2 



Y,n=o( X ( n ) - X ) COSUl(t n - t) 



Yln=0 COs2 - T ) 



T, n =o( x ( n ) ~ x ) sinw(i„ - r) 



and t is defined by the equation: 
tan(2wr) = 



J2n=o cos2w£„ 



(3) 



(4) 
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Fig. 1. A0224+671. Results for the 22 GHz (left panel) and 37 (right panel), daily averaged datasets. Sinusoids with 
a period equal to those provided by STIMA and listed in Tables 1 and 2 are overplotted as a visual aid to recognize 
periodicities. It needs to be stressed that these curves must be regarded as visual aids only: their amplitude is not by 
any means related to the signal. 




Fig. 2. 0945 + 408. Results for the 22 (left) and 37 GHz (right) M dataset. Sinusoids as in Fig. □ 



erties of the signal autocorrelation matrix (a.m.) which is 
the sum of the signal and the noise matrices. The p prin- 
cipal eigenvectors of the signal matrix allow the estimate 
of frequencies; the p principal eigenvectors of the signal 
matrix are the same of the total matrix. In order to ex- 
tract the principal components we use a robust nonlinear 
PCA Neural Network and we apply a modified version of 
MUSIC that can be applied directly on unevenly sampled 
data, to obtain the periodicities (cf. T99 and references 
therein). 

The STIMA process for periodicity analysis can be di- 
vided in the following steps: 

— Preprocessing: we first calculate and subtract the aver- 
age pattern to obtain a zero mean process (Karhunen 
& Joutsensalo 1994, 1995). 

— Neural computing: the fundamental learning parame- 
ters are: i) the initial weight matrix; ii) the number of 
neurons, which is the number of principal eigenvectors 
that we need, and therefore is equal to twice the num- 
ber of signal periodicities (for real signals); iii) a, the 
nonlinear learning function parameter; iv) the learning 
rate /j,. 

We then initialize the weight matrix W assigning the 
classical small random values. Otherwise we can use the 
first patterns of the signal as the columns of the matrix. 



Experimental results show, however, that even though the 
latter technique speeds up the convergence of our neu- 
ral estimator (n.e.), it cannot be used with anomalously 
shaped signals, such as stratigraphic geological signals. We 
use a simple criterion to decide whether the neural net- 
work has reached or not convergence: we control when the 
Power Spectrum of a modified MUSIC estimator is greater 
than zero. This modified MUSIC estimator can be written 
as: 

where w(i) is the i — th weight vector after learning, 
and ej is the sinusoidal vector: 

ef =[e t /,e},...,ef- 1 ] H (6) 

and to, tz,-i, are the epochs of the first L samplings. 
We then have the following general algorithm: 

- STEP 1: initialize the weight vectors Wo(i) Vi = 1, ...,p 
with small random values, or with orthonormalized signal 
patterns. Then initialize the learning threshold e and the 
learning rate //. Then reset pattern counter k = 0. 

- STEP 2: input the k - th pattern 

x fc = [x(k),...,x(k + N+l)} 
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Fig. 3. 1226+023. Panels a and b: results for the 22 GHz and the 37 GHz daily averaged radio curves, respectively. 
Panels c, d and e: the same for the 4.8, 8 and 14.5 GHz data, respectively. Sinusoids with a period equal to those listed 
in the tables are overplotted as visual aid (as in Fig.l). 



where N is the number of input components. 

- STEP 3: calculate the output for each neuron y(j) = 
w T (j)x t , Vi = 1, 

- STEP 4: Vi = 1, ...,p modify the weights 

w fe+ i(i) = w fe (j) + fi k g(yk(i))e k (i) 

- STEP 5: convergence test. If 

PM M-Efl\e?MiW >0 



in correspondence to the component frequency. Estimates 
are therefore related to the highest peaks. 

We note that MUSIC works on individual points, not 
on the basis of time intervals; it therefore does not need 
any kind of interpolation to deal with gaps in the data. 
Also, tests carried out by the Authors of this method 
(Tagliaferri et al. 1999) have shown that low-frequency 
drifts of the baseline flux of sources do not affect the de- 
rived periodicities. 

In the following discussion we shall make use of both 
STIMA and the Lomb Periodogram. 



, then GO TO STEP 7. 

- STEP 6: k = k + 1. GO TO STEP 2. 

- STEP 7. End. 3. Analysis of the Metsahovi radio light curves 

The frequency estimator MUSIC takes as input the 

weight matrix columns after the learning phase. The es- We have analyzed the 22 and 37 GHz, both daily and 

timated signal frequencies are obtained as the peak loca- weekly averaged, light curves of blazars monitored with 

tions of the function previously introduced. When / is the thc Metsahovi radio telescope, updated to the end of 2001 . 

frequency of the i — th sinusoidal component, f = fi , In order to establish whether a periodicity is real or 

we have Pm — > oo. In practice, we have a peak near and not we adopted the following set of rules: 
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Fig. 4. 2200 + 420. Panels have the same meaning as in Fig. [21 



— the periodicity derived by both STIMA and Lomb's 
method must be shorter than 1/2 of the time coverage 
(in order to avoid aliasing); 

— For a given object, periodicity must be found in both 
the 22 and the 37 light curves and the periods must be 
equal given a percentage of error (in the following we 
use a percentage of 10%) 

Of the 157 sources in the sample by Terasranta et al. 
(1998), 80 could not be analyzed due to either a very 
short time coverage or a too sparse distribution of the 
data points. Evidence of a periodicity was found for 5 of 
the remaining 77 objects (see Figs. ^Ell- The results for 
these, based on daily averaged fluxes, are given in Tablc^ 
where the left-hand side refers to 22 GHz data, the right- 
hand side to 37 GHz. The periodicity, given in units of 
10 3 days, is detected by both STIMA and Lomb's meth- 
ods and both methods yield, for each source, values very 
close to each other; the maximum difference is of ~ 8%. 
In the following, unless explicitly noted, we shall always 
use the STIMA estimate, which is less affected by spu- 
rious signals. If weekly, rather than daily, estimates are 
used, to increase the S/N ratio, the estimated periods are 



essentially unaffected (differences are at percent levels), 
even though the significance of periodicity is somewhat 
increased. 

Table [21 summarizes the results of the analysis of the 
UMRAO monitoring data for the five sources with de- 
tected periodicity, while a compendium of STIMA periods 
at the five UMRAO plus Metsahovi frequencies is pre- 
sented in Table [3 

For the sources 0224 + 671 (Fig. [TJ and 0945 + 408 
(Fig. 0) the number of data points is quite limited and, 
correspondingly, the periods are ill-defined. In the case of 
A0224 + 671 the analysis of UMRAO data yields signif- 
icantly different periods at 4.8 GHz compared to 8 and 
14.5 GHz; all of them are substantially higher (by roughly 
a factor of 2) than those found from the Metsahovi light 
curves. For 0945 + 408 (Fig.0) no periodicity was found at 
4.8 and 8.0 GHz, while an ill-defined period of 2690 days 
(STIMA), about a factor of 2 larger than found at the 
Metsahovi frequencies, could be present in the 14.5 GHz 
signal. These sources are good candidates for longer and 
more frequently sampled monitoring campaigns, particu- 
larly at the Metsahovi frequencies where the time span of 
light curves is relatively short. 
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Fig. 5. 2251 + 158. Panels have the same meaning as in Fig- EJ 




Fig. 6. OJ287. Panels a), b), c), and d): light curves at 22, 37, 14.5, and 4.8 GHz, respectively. The dot-dashed 
sinusoids have periods of about 20.9 yr (22 GHz), 22.3 yr (37 GHz), 20.6 yr (14.5 GHz), 25 yr (4.8 GHz). 
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Table 4. The U sample with source classifications by Do- 
nato et al. (2001) and Terasranta et al. (1998). In the lat- 
ter case BLO stands for BL Lac Object. Redshifts where 
taken from the SIMBAD astronomical database. 



Name 


Other name 


Donato 


Terasranta 


z 


0048-097 


OB -080 


LBL 


BLO 


0.2 


0133+476 t 


OC 457 


FSRQ 


HPQ 


0.859 


0202+149 t 


4C 15.05 


FSRQ 


HPQ 


0.405 


0235+164 t 


OD 160 


LBL 


BLO 


0.94 


0420-014+ 


OF -135 


FSRQ 


HPQ 


0.915 


0430+052+ 


3C 120 


FSRQ 


LPQ* 


0.0331 


0605-085 


OH -010 




HPQ 


0.872 


0607-157+ 


OH 112 






0.324 


0735+178+ 


OI 158 


LBL 


BLO 


0.424 


0814+425+ 


OJ 425 


LBL 


BLO 


0.2453 


0851+202+ 


OJ 287 


LBL 


BLO 


0.306 


0923+392+ 


4C 39.25 


FSRQ 


LPQ 


0.6948 


0945+408 


4C 40.24 


FSRQ 


LPQ 


1.252 


0954+658 




LBL 


BLO 


0.367 


1034-293 


OL -259 


FSRQ 


HPQ* 


0.312 


1055+018+ 


OL 093 


FSRQ 


HPQ 


0.888 


1127-145+ 


OM -146 






1.187 


1226+023+ 


3C 273 


FSRQ 


LPQ 


0.158 


1253-055+ 


3C 279 


FSRQ 


HPQ 


0.538 


1308+326+ 


OP 313 


LBL 


BLO 


0.996 


1418+546+ 


OQ 530 


LBL 


BLO 


0.151 


1510-089+ 


OR -017 


FSRQ 


HPQ 


0.359 


1538+149 


OR 165 


LBL 


BLO 


0.605 


1637+574+ 


OS 562 




LPQ 


0.7506 


1641+399+ 


3C345 


FSRQ 


HPQ 


0.594 


1642+690 


4C 69.21 




HPQ 


0.751 


1749+096+ 


OT 081 


LBL 


BLO 


0.322 


1749+701 




LBL 


BLO 


0.77 


1823+568 


4C 56.27 


LBL 


BLO 


0.6635 


1928+738 


4C 73.18 


FSRQ 


LPQ 


0.360 


2131-021 


4C -02.18 


LBL 


BLO 


1.285 


2134+004+ 


OX 057 


FSRQ 


LPQ 


1.931997 


2145+067+ 


OX 076 




LPQ 


0.99 


2155-152+ 


OX -192 


FSRQ 




0.672 


2200+420+ 


BL Lac 


LBL 


BLO 


0.0688 


2201+315+ 


4C 31.63 




LPQ 


0.298 


2223-052+ 


3C 446 


FSRQ 


HPQ 


1.404 


2251+158+ 


3C 454.3 


FSRQ 


HPQ 


0.859 


2254+074+ 


OY 091 


LBL 


BLO 


0.19 



+ Also in the Metsahovi database. 

+ Classified as FSRQ by Stickel et al. (1994). 

* Classified by Ghisellini et al. (1993). 



In the case of 1226 + 023 (Fig. |3J), the analysis of the 
UMRAO database yields well defined periods at all fre- 
quencies (4.8, 8, and 14.5 GHz). The periodic behavior is 
most evident at 4.8 and 8 GHz, but the estimated period 
at 8 GHz is lower by ~ 20% than at 4.8 GHz. At higher fre- 
quencies (14.5, 22 and 37 GHz) the periodic signal appears 
to be contaminated by an additional impulsive component 
which becomes dominant at 14.5 GHz. The estimated pe- 



riods at all frequencies, except 8 GHz, are consistent with 
each other to better than 10%. 

Periodicity is detected, by both the STIMA and the 
Lomb's methods, in light curves of 2200 + 420 at all fre- 
quencies (Fig. EJ), and period estimates are all consistent 
with each other , and close to the value obtained in the 
optical: 7.8 ± 0.2 yr, i.e. 2849 ± 73 days (Marchenko et al. 
1996; Hagen- Thorn et al. 1997). The periodicity is again 
most evident at 4.8 and 8 GHz, while at 22 and 37 GHz the 
residuals show the presence of an additional impulsive sig- 
nal. At 14.5 GHz, the evidence of periodicity is marginal, 
although the estimated period is very close to that found 
in the other bands. 

As shown by Fig. 2251 + 158 is probably our best 
case. Strong and consistent periodic signals are detected 
at all frequencies, most clearly at 22, 4.8 and 8.0 GHz. At 
all frequencies, however, an inspection of the the residu- 
als obtained by subtracting the first harmonic from the 
observed data, evidences an additional, impulsive, source 
of luminosity variations, over-imposed on the underlying 
periodic component. 

It may look surprising that our list of sources with ev- 
idence for periodicity does not include the most famous 
case, OJ287. The reason is that the estimated period does 
not meet our criterion of being shorter than 1/2 of the 
time coverage. In fact, we detect a possible modulation 
on a period of about 21-22 years atl4.5, 22, and 37 GHz, 
and on a longer period (about 25 years) at 4.8 GHz (see 
Fig. |SJ. These possible periodicities are about a factor of 
2 larger than found by analyses of optical data (Sillanpaa 
et al. 1988, 1996; Pietila, et al. 1999). An analysis of the 
residuals (after having subtracted the sinusoidal modula- 
tion) detects a significant periodicity of ~ 1.6 yr at 22, 37, 
and 14.5, in close agreement with that found by Aller et 
al. (1992) using a Scargle periodogram and by Hughes et 
al. (1998) using a wavelet transform. Data at 8 and 4.8 
GHz indicate (in the residuals) somewhat longer periods 
(1.7 and 2.05 yr, respectively). 

4. Phenomenology of variability properties 

The use of unbiased samples is essential for meaningful 
statistical analyses. The UMRAO group have monitored 
and analyzed the light curves of two complete samples, 
the Pearson-Readhead sample (Aller et al. 2003) and a 
BL Lac sample (Aller et al. 1999). 

The sample used in the present analysis is drawn 
from the complete catalog of sources brighter than 1 Jy 
at 5 GHz (Kiihr et al. 1981). We have have selected 
the Kiihr sources with UMRAO monitoring at 8 GHz 
for a total of > 3000 days, with gaps not exceeding 
200 days, all times being computed in the source frame 
(Source = t b server /(l + z)). Also, we have confined our- 
selves to sources classified either as Low-energy peak BL 
Lacs (LBLs, Padovani & Giommi 1995) or Flat-Spectrum 
Radio Quasars (FSRQ), thus excluding the lonely High- 
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Type 


N 


37 GHz 


22 GHz 


14.5 GHz 


8 GHz 


4.8 GHz 


LBL 


9 


0.74±0.27 


0.73±0.26 


0.69±0.19 


0.68±0.18 


0.59±0.16 


FSRQ 


16 


0.61±0.16 


0.55±0.15 


0.54±0.12 


0.50±0.11 


0.38±0.08 


HPQ 


9 


0.65±0.23 


0.60±0.22 


0.56±0.17 


0.51±0.16 


0.38±0.12 


LPQ 


7 


0.56±0.23 


0.48±0.21 


0.49±0.18 


0.45±0.17 


0.33±0.13 



Table 6. KS test on the distributions of variability indices. 





37 GHz 


22 GHz 


14.5 GHz 


8 GHz 


4i 


1 GHz 




d prob 


d prob 


d prob 


d prob 


d 


prob 


LBL vs FSRQ 


0.65 7.3 10 _a 


0.67 5.710 _a 


0.52 7.210" a 


0.61 1.4 10" 2 


0.67 


1.8 10~ 4 


HPQ vs LPQ 


0.32 0.73 


0.46 0.28 


0.39 0.33 


0.32 0.73 


0.28 


0.75 



Table 7. Mean values of structure function slopes. 



Type 


37 GHz 


22 GHz 


14.5 GHz 


8 GHz 


4.8 GHz 


all 

LBL 

FSRQ 


1.14 ±0.23 

0.98±0.37 

1.08±0.29 


1.18 ± 0.25 

0.95±0.37 

1.16±0.31 


1.17 ± 0.19 

1.05±0.29 

1.15±0.25 


1.17 ±0.20 

1.00±0.28 

1.19±0.26 


1.21 ±0.20 

1.05±0.29 

1.23±0.27 


HPQ 
LPQ 


1.05±0.40 
0.96±0.43 


1.09±0.41 
1.07±0.48 


1.08±0.35 
1.15±0.44 


1.18±0.38 
1.17±0.46 


1.14±0.37 
1.24±0.47 



energy peak BL Lac (HBL), namely Mkn 501, as well 
as the few sources classified as radio-galaxies or Seyfert 
galaxies. Whenever available, we have adopted the classi- 
fication given by Donato et al. (2001), supplemented with 
those by Terasranta et al. (1998), Ghisellini et al. (1993), 
and by Stickel et al. (1994). 

The final sample comprises 39 sources (24 FSRQs and 
15 LBLs), listed in Tabled (U sample). Of these sources, 
25 (9 LBLs and 16 FSRQs) were also extensively moni- 
tored by the Metsahovi group (UnM sample) . Of the 24 
FSRQs, 21 are classified as either HPQ (12) or LPQ (9); 
the acronyms stand for high- or low-polarization quasars, 
respectively. Redshift determinations are available for all 
sources. To characterize the variability properties at the 
various frequencies we have computed the variability in- 
dex (V.I.), the structure function (Simonetti et al. 1985; 
Hughes et al. 1992), and the distribution of intensity vari- 
ations. The variability index has been computed allowing 
for measurement errors following Aller et al. (2003): 

yj _ {S max - Q~S„ lax ) - [Smin + <J S min ) ^ 
{Smax — a S max ) + {Smin + 0"S mi „) 

Again following Aller et al. (2003) we have ex- 
cluded anomalously noisy data, i.e. data with as > 
max(0.1Jy, 0.035). 

4-1- Variability index 

Table [5] gives the mean values of the variability index 
(V.I.) and their errors at 4.8, 8, 14.5, 22, and 37 GHz for 
the 25 blazars (9 LBLs and 16 FSRQs) in the UnM sam- 
ple. Of the 16 FSRQs, 9 are HPQs and 7 are LPQs; we have 



also computed the V.I. for each of these sub-classes. At all 
frequencies we get systematically higher V.I.s for LBLs 
than for FSRQs (see also Fig. 0) , consistent with the re- 
sults by Aller et al. (1999). The Kolmogorov-Smirnov (KS) 
test shows that the differences among the two populations 
are significant: the probability that the two distributions 
come from the same parent population is ~ 3.1 x 10~ 3 , 
~ 2.4 x 10~ 3 , ~ 4.2 x 10~ 3 , ~ 8.6 x 10" 3 , and 9.1 x 10~ 5 , 
at 37, 22, 14.5, 8, and 4.8 GHz respectively (see Table®. 



The mean V.I. of FSRQ is higher at higher frequen- 
cies: it is ~ 0.38 at 4.8 GHz and becomes ~ 0.61 at 37 
GHz. Again, the KS test confirm that the effect is sta- 
tistically significant: the probability that the distributions 
at 4.8 and 37 GHz come from the same parent popula- 
tion is ~ 4.9 x 10~ 5 . The frequency dependence is less 
clear in the case of LBLs, which have large V.I.s at all 
frequencies. Still, their mean V.I. increases from ~ 0.59 
at 4.8 GHz to ~ 0.74 at 37 GHz; the significance of the 
difference between the distributions at the two frequen- 
cies is ~ 4.4 x 10~ 3 . The frequency dependence of the 
V.I. is consistent with the notion that, at lower and lower 
frequencies, an increasing fraction of the observed emis- 
sion is produced on larger and larger scales and is there- 
fore less variable on the timescales covered by monitor- 
ing programs. At all frequencies the mean V.I. of High- 
Polarization Quasars (HPQs) is slightly higher than that 
of Low-Polarization Quasars (LPQs); the difference, how- 
ever, has always a low statistical significance. 
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Table 8. Mean values of the log of the structure function timescales. 



Type 


37 GHz 


22 GHz 


14.5 GHz 


8 GHz 


4.8 GHz 


all 

LBL 

FSRQ 


0.04 ±0.07 
-0.03±0.07 
0.09±0.10 


0.10 ±0.08 

0.20±0.16 

0.07±0.10 


0.26 ±0.09 

0.20±0.13 

0.33±0.11 


0.27 ±0.09 

0.04±0.12 

0.41±0.13 


0.40 ± 10 

0.29±0.16 

0.46±0.13 


HPQ 
LPQ 


0.06±0.13 
0.17±0.17 


0.05±0.13 
0.13±0.15 


0.13±0.10 
0.46±0.22 


0.08±0.12 
0.61±0.26 


0.19±0.14 
0.64±0.26 




Fig. 7. Variability index distributions at different frequencies for the UflM sample. 



4-2. Structure function 

Another useful method for quantitatively investigating 
variability properties is the so called structure function 
(S.F.) analysis (Simonetti et al. 1985; Hughes et al. 1992). 
We have computed the quantities commonly used to char- 
acterize the first order S.F., defined as D(t) — {[S(t) — 
S(t+r)] 2 ), namely its slope (<i log Z?/<i log r) and the time- 



scale, i.e. the time lag of the turnover of the S.F., which 
may measure the minimum time scale of uncorrelated be- 
havior or, in the case of flicker noise, the minimum time 
scale in the distribution of response times (Hughes et al. 
1992). 

In Fig.[S]we show the structure functions at 22 GHz for 
the 4 sources for which we have the best evidence of peri- 
odicity. The value of the derived period is indicated by the 
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Fig. 8. Structure functions at 22 GHz of the 4 sources for which we have the best evidences of periodicity. The values 
of the periods are indicated by the vertical lines, while the dotted horizontal lines show the variance of the process. 



vertical line, while the dotted horizontal line corresponds 
to the variance of the process. 

The distributions of slopes and of time-scales (in the 
source frame) are shown in Figs. El and respectively. 
Panels corresponding to frequencies of up to 14.5 GHz 
refer to the 39 sources in the U sample, those at higher 
frequencies to the 25 sources in the UflM sample. Mean 
values and dispersions are given in Tables and [SJ Slopes 
and time-scales (or lower limits) could be determined for 
essentially all sources. The exceptions (sources with very 
irregular S.F.) are just one at 14.5 and 8 GHz, and 3 at 4.8 
GHz. For the sources (2 at 37 and 22 GHz, 10 at 14.5 GHz, 
11 at 8 and 4.8 GHz) not showing a plateau we adopted 
the time base of the data as a lower limit to the time-scale. 
We further have a small number of sources with a plateau 



(1 at 37 GHz, 4 at 14.5 GHz, 1 at 8 GHz, and 2 at 4.8 
GHz) or a change of slope (1 at 37 GHz, 2 at 22 GHz, 5 at 
14.5 GHz, and 4 at 8 and at 4.8 GHz) at intermediate time 
lag, which may signal the transition between two different 
processes; in these cases we have adopted the longer time- 
scale. For sources with a change of slope we have chosen 
the steeper value. 

The average time-scale tend to be somewhat longer 
(although the statistical significance of the difference be- 
tween 37 and 4.8 GHz is, for the full sample, only 1.2 x 
10~ 2 ), and the fraction of sources with time-scale exceed- 
ing the data time span is larger at lower frequencies. The 
indication, reported by Hughes et al. (1992), that LBLs 
have, on average, shorter timescales than FSRQs is not 
statistically significant in the source frames (see Table |HJ). 
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0.5 1 1.5 2 0.5 1 1.5 2 

SF Slope SF Slope 



Fig. 9. Distributions of structure function slopes at different frequencies. The results at frequencies < 14.5 GHz refer 
to the full U sample, while at higher frequencies we are forcefully limited to the 25 sources also monitored at Metsahovi 
(UHM sample). 



There is a marginal indication (significance of 2.4 x 1CP 2 
based on the KS test at 8 GHz) of longer timescales for 
LPQs, compared to HPQs. 

No statistically significant variations with frequency of 
the slope distribution are detected. A visual inspection of 
Fig. indicates that FSRQs tend to have steeper slopes 
than LBLs. The average slope of LBLs is ~ 1, the value 
corresponding to shot noise. In the the case of FSRQs, the 
distributions are somewhat broader and extend to values 
> 1.5, particularly at lower frequencies, consistent with 
the previous findings by Hughes et al. (1992). The KS 
test yields a probability of 1.2 x 10~ 2 at 4.8 GHz and of 
4.3 x 10 -2 at 8 GHz that the distributions for the two 
populations come from the same parent distribution. No 



statistically significant differences are found at higher fre- 
quencies. The distributions of LPQs and HPQs are con- 
sistent with the same parent population at all frequencies 

Possible hints of a correlation of time-scales with slopes 
were noted by Hughes et al. (1992). A Pearson's test how- 
ever does not detect any significant correlation among 
these quantities either for the entire sample or for any 
sub-population. 

4-3. Distribution of intensity variations 

As shown by Fig. El the distribution of log(5/ S), where S 
is the instantaneous flux density at a given frequency and 



A. Ciaramella et al.: A multifrequency analysis of radio variability of blazars 



13 



6 

4 

2 


6 

4 

2 


8 

6 
z 4 
2 

10 
8 
6 
4 
2 

e 

6 
4 
2 





22 GHz 



~| I I I I | I 

37 GHz 

LBL 

_ _ FSRQ r 



rLnJ 



i i 



14.5 GHz 



8 GHz 




i - 



i j 

I E 



L4.8 GHz 




-0.5 







0.5 



log(Time Scale, yr) 



6 r 

4 r 

2 r 

r 

6 r 

4 r 

2 r 

- 

8 - 



HPQ 
LPQ 



I - 



i ' I ' 




-0.5 







0.5 



log(Time Scale, yr) 



Fig. 10. Distributions of structure function timescales at different frequencies. The samples used are the same as in 
Fig. El 



S is its average value, broadens with increasing frequency, 
consistent with the notion the the variability amplitude in- 
creases with frequency (Impey & Neugebauer 1988). The 
frequency dependence of the variance, tr 2 , of the distri- 
bution is approximately linear. A fit, shown as a dotted 
line in the upper part of the panel on the lower right-hand 
corner of Fig. 1111 is given by: 

a 1 ~ 1.846 x KT 4 i/ G Hz + 0.01837 . (8) 

The shape of the distribution also changes somewhat with 
frequency. At high frequencies (22 and 37 GHz), the dis- 
tribution of intensity fluctuations approaches a log-normal 
distribution, while at lower frequencies, the distribution of 
log(Sy S) is better described by a Laplace or by a Cauchy 
distribution. The frequency dependencies of the skewness 



and of the kurtosis are shown in the middle and on the 
bottom part of the panel on the lower right-hand corner 
of Fig. El 

Variability enhances the bright portion of the luminos- 
ity function and of source counts. We have estimated its 
effect on the 37 and 100 GHz counts of FSRQs and of 
BL Lacs by convolving the epoch-dependent luminosity 
functions given by Maraschi & Rovetti (1994) and Urry & 
Padovani (1995), respectively, with the distribution func- 
tion of intensity fluctuation shown in the corresponding 
panel of Fig. ^2 at 37 GHz, and with a Gaussian with 
variance given by Eq. JSJ at 100 GHz. As illustrated by 
Fig. 1121 the enhancement of the counts is of about 20- 
30% at the brightest flux densities. It is thus very likely 
that a substantial fraction of blazars picked up by the 
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Fig. 11. Probability distribution functions of \og(S/S at different frequencies and their Gaussian fits. The panel in 
the lower right-hand corner shows, as a function of frequency, the variance (multiplied by 100), the skewness (xlO), 
and the kurtosis. 



high-frequency, shallow surveys carried out by the WMAP 
(Bennett et al. 2003) and Planck satellites 1 are in a flar- 
ing stage. 



5. Discussion and conclusions 

Combining the long-term monitoring databases of the 
University of Michigan Radio Astronomy Observatory 
(UMRAO) and of the Metsahovi Radio Observatory it 
is possible to investigate the radio variability properties 
of a rather rich sample of FSRQs and LBLs over about a 
decade in frequency. 



1 astro.estec.esa.nl/SA-general/Projects/Planck/ 



Our approach has followed different lines of investiga- 
tion. First, we have addressed the still controversial issue 
of the existence of periodicities in the radio light curves. 
In addition to the well known Periodogram method, mod- 
ified to take into account the effect of uneven sampling, we 
have applied a more effective technique, never used before 
in this context, based on the analysis of the signal auto- 
correlation matrix. We have analyzed the Metsahovi light 
curves looking for periods shorter than 50% of the maxi- 
mum admissible period (to avoid aliasing) and detectable 
(with the same value within the errors) at both 22 and 37 
GHz. We have found evidences, confirmed by both tech- 
niques, for periodicities in the light curves of 5 (0224+671, 
0945 + 408 1226 + 023, 2200 + 420, and 2251 + 158) out 
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Fig. 12. Effects of variability on counts of FSRQs (solid 
lines) and of LBLs (dashed lines). For each population, the 
upper curve refers to 37 GHz, the lower one to 100 GHz. 
We have adopted a Gaussian distribution for log(Sy S with 
variance given by Eq. JHJ. Shown are the percentage differ- 
ences between the counts estimated allowing for variability 
and those assuming that each source keeps at its average 
flux S. 



of the 77 Metsahovi sources having sufficient time cover- 
age for a meaningful analysis to be carried out. For the 
last three of these sources, consistent periods are found 
also at the three UMRAO frequencies, while in the case 
of 0224 + 671 and 0945 + 408 (which were monitored for 
a relatively short time at Metsahovi) UMRAO data give 
hints of ill-defined periods a factor of 2 larger than at the 
Metsahovi frequencies. 

The Lomb periodogram method allows us to test quan- 
titatively the significance of the detected periodicities. As 
shown by Scargle (1982), the false-alarm probability, i.e. 
the probability that, if we scan M independent frequen- 
cies, none has spectral power normalized to the variance 
of the data larger than z, under the null hypothesis that 
the data are independent random Gaussian values, is: 

P(> z) = 1 - (1 - e- z ) M . (9) 

An accurate evaluation of P(> z) is complicated by the 
difficulty of estimating M (Press et al. 1996, p. 570). How- 
ever, for any reasonable choice of M , P(> z) turns out 
to be extremely small (and therefore the significance of 
the peak in the power spectrum corresponding to the de- 
tected period is very high) for 1226 + 023, 2200 + 420, and 
2251 + 158 (z amounts to several tens at all frequencies). 
On the other hand, the false alarms probability may be 
quite significant for 0224 + 671 and 0945 + 408. As a fur- 
ther test, we have applied the Kolmogorov-Smirnov test 
to examine the consistency of the data distribution with a 



Gaussian random process. This possibility is ruled out for 
1226 + 023 and 2251 + 158, but not for the other sources. 

We have also investigated the variability index, the 
structure function, and the distribution of intensity varia- 
tions of the most extensively monitored sources. We con- 
sidered UMRAO sources classified either as Low-energy 
peak BL Lacs or Flat-Spectrum Radio Quasars, moni- 
tored at 8 GHz for a total of > 3000 days, with gaps 
not exceeding 200 days, all times being computed in the 
source frame. The sample comprises 39 sources (24 FS- 
RQs and 15 LBLs), 25 of which (9 LBLs and 16 FSRQs) 
were also extensively monitored by the Metsahovi group. 
We have found a statistically significant difference in the 
distribution of the variability index for LBLs compared to 
FSRQs, in the sense that the former are more variable. 
This difference may help shedding light on the relation- 
ship between the two blazar sub-classes and is consistent 
with LBLs having a higher ratio of beamed to unbeamed 
emission than FSRQs (Fan 2003). 

For both populations the variability index steadily in- 
creases with increasing frequency. The distribution of in- 
tensity variations also broadens with increasing frequency, 
and approaches a log-normal shape at the highest frequen- 
cies. Variability enhances by 20-30%, at bright flux den- 
sities, the high frequency counts of extragalactic radio- 
sources, such as those which will be carried out by the 
WMAP and Planck satellites. 

Since the sample is not complete, one should worry 
about possible selection biases. To test this possibility we 
have extracted from it a sub-sample 80% complete to a 
flux limit of 1.2 Jy over the three areas 16 h < a < 20 h , 
53° < 5 < 75° ; 21 h < a < 23 h , -8° < 5 < 45° ; 
8 h < a < U h , 30° < 5 < 45° . Such sub-sample com- 
prises only 16 objects (6 LBLs and 10 FSRQs) and is 
thus too small for our purposes. On the other hand, we 
did not detect significant differences in the mean prop- 
erties of sources in the sample of 39 compared to those 
in the almost complete sub-sample. Also the trends with 
frequency and the differences among sub-populations are 
still present. We are therefore confident that, although the 
sample is not complete, it is essentially unbiased. 
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